Collapsing Bose-Einstein condensates beyond the Gross-Pitaevskii approximation 
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We analyse quantum field models of the bosenova experiment, in which **^Rb Bose-Einstein con- 
densates were made to collapse by switching their atomic interactions from repulsive to attractive. 
Specifically, we couple the lowest order quantum field correlation functions to the Gross-Pitaevskii 
function, and solve the resulting dynamical system numerically. Comparing the computed collapse 
times with the experimental measurements, we find that the calculated times are much larger than 
the measured values. The addition of quantum field corrections does not noticeably improve the 
agreement compared to a pure Gross-Pitaevskii theory. 
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I. INTRODUCTION 

Dilute gas Bose-Einstein condensates (BECs) display 
not only interesting classical nonlinear field physics, but 
also quantum field physics, such as the Mott insulator 
transition These BECs are especially interesting be- 
cause the weak atomic interactions make the physics suf- 
ficiently simple that we expect to get quantitative agree- 
ment between experiments and theory. In this paper 
we analyse one of the more straightforward aspects of 
the JILA bosenova experiment 0: the time to collapse. 
It has previously been shown that the Gross-Pitaevskii 
(GP) theory substantially overestimates these collapse 
times • We show that adding the lowest order quantum 
field corrections does not improve the situation. Thus an 
open question remains: what explains the discrepancy 
between theory and experiment? 

The collapse experiments at JILA 0, 0] represent the 
first completely controlled studies of unstable BECs. 
Earlier experiments on unstable ^Li condensates |^ Q 
could not provide such a well defined initial state for the 
collapse since the condensate was continuously fed by a 
thermal cloud and collapse was therefore triggered in a 
random fashion. Nonetheless these were the first exper- 
iments to verify theoretical predictions about the criti- 
cal atom number in attractive BECs. In the experiment 
described in Ref. 0, a stable ^^Rb condensate was pre- 
pared with vanishing interactions, and then the scatter- 
ing length was switched to a negative (attractive) value 
using a Feshbach resonance. 

The experiment demonstrated a number of compli- 
cated and interesting features of the condensate collapse, 
which was controlled via the use of a Feshbach reso- 
nance. We have extended exis ting theoretical studies 

saiiiiaiiiiiiiiiiiiimiiiiiiiaiii by 

applying the dynamical Hartree-Fock Bogoliubov (HFB) 
equations with realistic experimental parameters. 
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This paper is organized as follows: In the first part 
we will give an overview of the JILA experiment and the 
associated theoretical literature. We then introduce our 
quantum field models, and finally present the results of 
our simulations. 



A. Experiment 

The collapsing condensate was observed to lose atoms 
until the atom number reduced to about the critical value 
below which a stable condensate can exist j^- The de- 
pendence of the number of remaining atoms on time since 
initiation of the collapse TcvoIvc was measured for the case 
of an initial state with A^init ~ 16000 atoms and repul- 
sive interaction corresponding to Cinit = +7ao, where ao 
is the hydrogen Bohr radius. The onset of number loss 
is quite sudden, with milliseconds of very little loss fol- 
lowed by a rapid decay of condensate population (within 
0.5 ms) after which the condensate stabilizes again. This 
behavior results from the scaling of the loss rate with 
the cube of the density, the peak value of which rises as 

1/ (^collapse 

- t) near the collapse point . This allows 
a precise definition of the collapse time tcoUapse, the time 
after initiation of the collapse up to which only negligible 
numbers of atoms are lost from the condensate. Another 
quantitative result of the experiment is the dependence 
of icoiiapsc on the magnitude of the attractive interaction 
that causes the collapse, parametrised by the (negative) 
scattering length flcoUapsc- These measurements are per- 
formed from an initial state with A'init = 6000 atoms in 
an ideal gas state (with interaction between them tuned 
to zero) . The icoiiapsc datapoints presented in the original 
paper have undergone one revision of their acoUapso values 
by a factor of 1.166(8) due to a more precisely determined 
background scattering length Although the main fo- 
cus of this paper shall be on the collapse time, we mention 
two other striking features of the experiment: the appear- 
ance of 'bursts' and 'jets'. One fraction of the atoms that 
are lost during the collapse is expelled from the conden- 
sate at quite high energies (~100 nK to ~400 nK, while 
the condensate temperature is 3 nK); this phenomenon 







Ninit 




case 


(i) 


16000 


+7ao 


case 


(ii) 


6000 





case 


(iii) 


4000 


+2500ao 



TABLE I: Initial states for different collapse scenarios. 

was referred to as 'bursts'. Finally, when the collapse was 
interrupted during the period of number loss by a sud- 
den jump in the scattering length, another atom ejection 
mechanism was observed: 'jets' of atoms emerge, almost 
purely in the radial direction and with temperatures a lot 
lower than that of the bursts (a few nK). In addition to 
the two initial scenarios described above there have been 
unpublished measurements on a third initial state with 
Nin^t = 4000 and ai„it = -l-2500ao HH. To simphfy refer- 
ence to the three different initial states for the collapse, 
we label them as in Tab. |l| 

B. Theory 

1. Collapse time 

The unstable attractive BEC and its loss of atoms dur- 
ing collapse can be modeled numerically with the Gross- 
Pitaevskii equation (GPE) if a phenomenological three 
body loss term is included. Numerical solutions to this 
equation for exact experimental parameters and geom- 
etry have been reported by Adhikari 0, 0; Saito and 
Ueda H mill Santos and Shlyapnikov O, Sav- 
age, Robins and Hope Q as well as Bao, Jaksch and 
Markovich ■ Where they consider the initial situation 
of case (ii) , these solutions agree on collapse times which 
qualitatively describe the experiment but systematically 
exceed the experimental values. Note that while the rem- 
nant atom number after collapse depends strongly on the 
three-body recombination rate K^, which is not well de- 
termined near the Feshbach resonance, and for which dif- 
ferent values were employed in the references above, the 
collapse time in case (ii) does not vary much for exper- 
imentally reasonable values of K3 (3|. Fig. ^ shows a 
comparison between GP theory (•) and the revised ex- 
perimental data (x) for case (ii). Also shown are the 
GP results for spherical symmetry, for comparison with 
our quantum field theory calculations. The change in 
Ocoiiapso by a factor of 1.166(8) due to the revision has 
been included and moved the experimental and theoret- 
ical points closer together but they still do not agree. 

In 0, 0, Q| this situation is accepted as agreement 
between theory and experiment. But if the GP equa- 
tion indeed contains all the essential physics to describe 
^collapse, one should expect better agreement since all cru- 
cial parameters of the model are experimentally fixed. In 
addition there exists another case where the GPE can ap- 
parently not account for ^collapse accurately, and that is 
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FIG. 1: Experimental and numerical results for the collapse 
time tcoiiapse versus scattering length flcoiiapse after a switch 
from a = to (Jcoiiapsc in case (ii). The experimental points 
( X ) and their error bars are from Fig. 2 of _2j taking into ac- 
count the revision of scattering length values in |2C| . Numer- 
ical results are given for exact, cylindrical geometry (•) and 
for spherical geometry (*) with trap frequency equal to the 
geometric mean of the experimental values. K3 = 1 x 10~^^ 



the above mentioned case (iii) [l3- One reason for the 
controversy in the literature might be the fact that the 
GPE describes collapse times and remnant numbers well 
for case (i) 0, 0, S Q • The condensates for cases (i) 
and (ii) are similar in peak density, ruling out a simple 
qualitative difference of the densities involved as an ex- 
planation for this distinct behaviour. 

Complementary to the numerical approach to the 
problem, one can attempt to find an analytic approxi- 
mation that describes the essential features of the BEC 
collapse. Metens et al. 0| employ a variational ap- 
proach assuming a Gaussian condensate shape. For val- 
ues ||a| — Iflcrl |/|acr| << 1, where Ccr is the critical scat- 
tering length above which the BEC does not collapse, 
they obtain icoiiapse ~ (|a/acr| — 1)~^^^, in agreement 
with But this parameter regime does not apply to 
the experiment, and it turns out that the scaling law 
becomes modified to tcoiiapse (|a/acr-| — 1)""'^^^ for scat- 
tering lengths away from the critical value . Calzetta 
et al. |0| derive a scaling law for the collapse time from 
the growth of unstable perturbations of the atom field. 
Their result agrees with the one cited above, which seems 
to capture the shape of the function correctly. However 
both methods are unable to account for the proportion- 
ality constant, as both groups fit this parameter or the 
initial state to obtain agreement with experimental data. 

In summary, mean field theory can describe the shape 
of the dependence of tcoiiapse on acoUapse correctly, but for 
case (ii) the real collapse seems significantly accelerated 
compared to these predictions. Naturally, one would sus- 
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pect deviations from GP theory to be the cause of this 
phenomenon IT^ . 

Yurovsky |22l | investigates quantum effects in homoge- 
neous BECs with attractive interaction analytically, us- 
ing an approximate solution to the operator equation for 
the condensate fluctuations. Results from regarding 
the growth of an initially seeded unstable momentum- 
mode have been experimentally verified 23]. Yurovsky 
also derives an expression for the time evolution of con- 
densate depletion. 



2. Jets 

Another prominent feature of the JILA experiment was 
the appearance of radial atomic jets when the collapse 
was suddenly interrupted. Numerical simulations were 
able to reproduce them using the GPE The 
phenomenon can be explained as interference of atom 
emission from highly localized density spikes that are 
formed during collapse _9i. The simulations seem to re- 
produce jet-shapes, anisotropy and dependence on ^evolve 
(Fig. 5 of 2]). By contrast Calzetta et al. present a 
mechanism for jet production that is essentially a quan- 
tum effect. Their model reproduces the features of the 
experimental data for atom numbers in jets (Fig. 6 of 0) 
for the range of validity of the approximations involved. 
Regarding this, the numerical simulations in 

Bill yield 

a similar (partial) agreement with experiment, but do not 
agree with each other. Another feature of the jets in the 
collapse experiment is the existence of large fluctuations 
in the atom number in jets. Calzetta et al. |l7j | attribute 
this to a large number uncertainty in squeezed quantum 
states, while Bao et al. reproduce these fluctuations 
by slight variations in the position of the initial state con- 
densate wavefunction. The existence of jets can thus be 
predicted by GP theory alone, but it is not yet clear what 
role quantum effects play. 



of Qcoiiapse, One caunot obtain simultaneous agreement 
for cases (i) and (ii) described in Ref. 0|- Furthermore, 
the remnant atom number after collapse iVj-cmn also de- 
pends on and the values for that are necessary 
to obtain high energy bursts in agreement with the ex- 
periment yield A^remn that disagree with the experiment. 
Calzetta et al. [T^ also suggest a quantum treatment 
of the bursts but do not put forward any qua ntitative 
predictions. Milstein, Menotti and Holland [13 present 
a model calculation that shows that atom-molecule cou- 
pling can yield sufficiently energetic atoms, however they 
do not examine the experimental parameter regime. In 
this paper, we will follow the formalism developed by Mil- 
stein et al. but place a lower emphasis on the molecular 
field. 

An alternate model that describes physics beyond the 
GPE has been proposed by Duine and Stoof [1^ [Tjl and 
applied to the collapse case (iii). They suggest elastic 
collisions between condensate atoms as the origin of the 
burst, and an additional channel by which the conden- 
sate can eject atoms. One of the collision partners gets 
ejected into the noncondensate in the process, while the 
other one is stimulated back into the condensate. The 
model yields results that describe well the experimental 
data for the case (iii) collapse. But case (ii) cannot be 
treated with a Gaussian variational ansatz 18], and it is 
not clear how to incorporate these elastic collisions into 
a treatment that is focussed on the backreaction of un- 
condensed atoms on the condensate. We will comment 
further on the relation to our model in the last section. 

A numerical study that accurately reproduces the 
bursts and all other results of the experiment therefore 
does not yet exist. 



II. THE HARTREE-FOCK BOGOLIUBOV 
MODEL 

A. Quantum corrections 



3. Bursts 

Finally, one has to explain the bursts. These are atoms 
that remain in a detectable, trapped state during collapse 
but are found performing large amplitude oscillations in 
the trap. Numerical solutions of the GPE generally show 
atoms being ejected from the condensate at the time of 
collapse. In order to be candidates for the burst atoms in 
the JILA experiment, they need to have rather high ener- 
gies ^ 100 nK (compare Fig. 4 of [2|). The ejected atoms 
within the condensate gain their energy due to a loss in 
(negative) mean field energy by three body recombina- 
tion, which is sensitive to the three body loss rate K3. 
Since K3 is experimentally not very well constrained in 
the parameter regime in question |24|, it is possible to 
reproduce the atom burst by matching on the en- 
ergy spectrum But in Ref. it is shown, that 
even if is fitted to the energy spectrum for each value 



To go beyond mean field theory we start from the usual 
Hamiltonian for a many body system that interacts by a 
contact potential of the form V{x — x') = C/q'^*-'^'' (x — x'). 
In the following, ^^(x, i) will denote the field operator 
that annihilates an atom at position x. 

H= ydx^'t(x)7?0a(x)^'a(x) 

+ dx*t(x)*t(x)|-,(x)*,(x), (1) 

where 

F(x) = \n{i^\x^ + uj\y^ + cjjjz^). (2) 
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Here m is the atomic mass (1.41 x 10~^^ kg for ^^Rb), 
Uo the atom-atom couphng constant and uj±, ui\\ de- 
note the trapping frequencies in the transverse and ax- 
ial directions. Now we derive the Heisenberg equa- 
tion for the field operator, and subsequently decom- 
pose \['q(x, i) into a condensate part (/)a(x, i) and quan- 
tum fluctuations x(x, i), such that '^a = 0a + X ^nd 
(^) — (j)a- We will describe the fluctuations in terms 
of their lowest order correlation functions: the normal 
density G'Ar(x, x') = (x^(x')x(x)) and anomalous den- 
sity Ga(x, x') — (x(x')x(x)). In deriving the dynamical 
equation for the condensate we will factor the expecta- 
tion values in accordance with Wick's theorem, therefore 
assuming that the system is in a Gaussian quantum state 
(i.e. a coherent state or even a squeezed state) and 
obtain: 



dt 



-^vi + F(x) + Uo\M^)V ] M^) 



+ 2UoM^)Gn (x, x) + UoC^I{x)Ga (x, x) 



(3) 

This modified Gross-Pitaevskii equation now contains 
the interaction of the uncondensed component with the 
mean field. We also model three-body loss from the con- 
densate phenomenologically, by including the following 
term in equation Eq. Ipj)l: 



-*-i^3|0a(x)|'*(/..(x). 



(4) 



As the next step we derive the evolution equation for 
Gn/a from the Heisenberg equations for the operators 
X^{x')x{x) and x{x')x{x) which gives us: 

{Haa{x) + Hoa{x'))GA (X,X') + 2Uo{ |0a(x)|' + \M^')f 

+ Gn (x,x) +Gn (x',x'))Ga (x,x') 

+ Uo{M^fG% (x, x') + M^TGn (x, x') 

+ Ga (x, x) (x, x') + Ga (x', x') Gn (x, x') ) 

+ Uo{M^? + Ga (x,x) )<5(3)(x - x'), (5) 

^h^^^Lp^ = {[xH-'m^lH]) = 

(i?Oa(x) - Hoa{x'))GN (x, x') + 2[/<,( |0,(x)|' - \M^')f 

+ Gn (x, x) - Gat (x', x') ) Gn (x, x') 

+ Uo {Ga (x, x) G*a (x, x') - G*a (x', x') Ga (x, x') ) 

+ Uo{M^)'G*A (x, x') - Ki^TGA (x, x') ) . (6) 

Equations l|3Il-® constitute the so called dynamical HFB 
equations and are identical to those of Ref. ,lSi] , without 



the molecular field (see below) . The calculation of quan- 
tum corrections in a local field theory naturally leads 
to divergences. Therefore a renormalization of the the- 
ory is necessary. As the renormalization condition we 
demand that the amplitude for low energy atom scat- 
tering is given hy U = 47r fi^a/m where a is the s-wave 
scattering length Ultraviolet (high momentum) di- 
vergences that arise when this amplitude is determined 
from the Hamiltonian are regularized by a momen- 
tum cut-off K. This yields a relation between the phys- 
ical interaction strength U and the parameter Uo in the 
Hamiltonian. Following Kokkelmans et al. |2^ : 



Uo 



u 



1-aU' 



mK 



(7) 



This process also provides a natural interpretation of the 
delta function in Eq. jSJ in terms of its truncated Fourier 
transform. Some subtle issues connected to the renormal- 
ization procedure will be discussed in the final section of 
this paper. The reason that quantum field effects might 
improve the agreement in Fig. ^ is the statistical factor 
of two in the interaction between the atomic mean field 
and the uncondensed component in Eq. (PI . If a sufficient 
number of atoms are transferred to uncondensed modes 
during collapse, it is expected to occur earlier. 



B. Resonance theory 

In the vicinity of a Feshbach resonance the atomic in- 
teractions are strongly energy dependent. To investigate 
effects of this energy dependence it is possible to model 
the Feshbach resonant interaction by including a molec- 
ular field (x) in the model • The additional terms 
in the Hamiltonian are then: 



/ 



dx *t„(x) + 2T/(x) - ) ^^(x) 

I J dx*L(x)^a(x)^a(x)+ff., 



(8) 



and the renormalization formalism presented in |27j| now 
connects the atom- molecule coupling g, molecular detun- 
ing v and background atom-atom scattering Uhg to the 
effective interaction Ucs that the atoms experience due 
to their coupling to the molecular field by: 



Uo^ 



1 - aUbg' 
9 

I - aUbq' 



9 



99o 

Vn = V — a = V — a- , 

2 2(1 -a [/bg)- 



9 



t/cff - C/bg I 1 



2C/bgJ^ 



(9) 
(10) 

(11) 
(12) 
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Note that the parameters Uhg,g,i^ are physical parame- 
ters, distinct from the bare parameters UoTgoi^^o^ which 
are just the parameters in the Hamiltonian and depend 
on the chosen cut-ofF. The atom-molecule coupling is 
related to the width of the Feshbach resonance k by: 

9 = V^^^bg = V^A/ibg , (13) 

where F is the width in terms of the magnetic field [2(tI |. 
and A/ibg is the difference in magnetic moments of the 
open and closed channels |2||. The coupling is therefore 
experimentally determined, as is the background scatter- 
ing length C/bg >2fiJ- 

In deriving the equations of motion for the cou- 
pled atom-molecule system quantum fluctuations of 
the molecular field are neglected, that is we assume 
(^',„(x)} = (?!)„i(x) which yields: 

+ 2UoGn (x, x) cpa{^) -f UoGa (x, x) 0:(x) 

+ 5o0™(x)(^:(x), (14) 

+ |^(G^(x,x)+</)„(x)2). (15) 

For the modifications of the remaining two Hartree-Fock 
Bogoliubov equations 1)516(1 by the coupling to molecules 
we refer to Ref. [iflj . 

III. RESULTS 

A. HFB Approach 

While the GPE can be tackled even in a completely 
asymmetric situation 's'l , doing so for the set of equations 
(inj-® even with cylindrical symmetry would present a 
serious numerical problem since the dimensionality of the 
correlation functions would only reduce from six to five. 
The resulting requirements for memory and computation 
time would exceed those for spherical geometry by two 
orders of magnitude. One problem is, that the violent 
nature of the collapse and the resulting spreading of the 
momentum space wavefunction do not allow radial grids 
smaller than 128 points, and our experience is that the 
minimum number of grid points to describe an angular 
variable lies around 8. Higher numbers are preferable. 
We are therefore forced to model the experiment with 
a spherically symmetric model. The trap frequency is 
chosen to equal the geometric mean of the experimen- 
tal trap frequencies: to = (a;^W||)^/'^. With wj_ =17.5 
Hz and uj\\ =6.8 Hz from 0, oj =12.8 Hz. To determine 
the influence of geometry on the collapse, we numerically 



solved the spherically symmetric GPE for this trap pa- 
rameter. The results are included in Fig.n^^nd show that 
the spherical situation gives quite a good estimate of the 
complete three dimensional GPE collapse time. In order 
to simplify equations ((SJ and © for the case of spherical 
symmetry we work in coordinates: 

i?=|x|, i?' = |x'|, /3 = cos6', (16) 

where is the angle between x and x'. Due to the restric- 
tion to spherical symmetry the correlation functions will 
not depend on the overall orientation of x and x', or on 
the azimuthal angle of x' around x, reducing the dimen- 
sionality of Ga and Gn from six to three. Following 
the (3 dependence is expanded in terms of Legendre poly- 
nomials Pn (/?) . Numerically, we model the delta function 
in Eq. (jSJ by a step function on one spatial grid point. 
The corresponding cut-off in Eq. {TJ is thus K — tt/AR, 
where AR is our grid spacing. Since we calculate spatial 
derivatives with the EFT all radial grids need to incor- 
porate an unphysical negative range. We now define the 
functions: 

(17) 

and expand the correlation functions as: 

M-l 

Gn {R, R', f3)=J2 (18) 

n=0 
M-l 

Ga{R,R'.P) = {R.R')Pn{P), (19) 

n=0 

where M is the number of Legendre Polynomials em- 
ployed. We use the additional notation: 

Gm/a{R) = g^m(R, R) = E g%a(R^ R)- (20) 

The final form of our equations reads: 

.^dUR) _ 



(-2;;^ai?^ + ^(^)-^2i?i^"(^)iT^(^) 

+ ^{2Gn (R) + \4>a{R)\') MR) 

+ ^Ga{R)K{R). (21) 
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for the atomic condensate. The expansion coefficients 
obey the following equations of motion: 



ih 



dt 



1 1 



V{R) - V{R') + 2Uo 



\^a{R)?+GN (R) 



- 2U, 



\^a{RT + GM {R') 



i?2 

{R,R') 



4>a{Rf +Ga {R) ^{n)* ,r. n/^ 

Uo ^2 y^iti) 



(22) 



ih 



dG^2^ {R, R') _ 



V{R) + V{R') + 2U, 

\4>aiR')\' + GN {R') 



\MR)? + Gn (R) 



i?2 



2Uo^ 



i?'2 



G^A^ {R,R') 



Uo 



MR?+Ga (R) Mn)* 



i?2 



( n ril\ 



Uo 



4>a{R'?+GA {R')Mn) 



U 



i?'2 

o[MRf + GA (R) 



G^''' {R, R') 



2n + l 



S{R-R'). 



(23) 



47ri?2 

As initial checks of our code we reproduced the results 
of Milstein et al. and Holland et al. using the 

molecular field and the equations of Ref. |19| |. We then 
computed a solution to Eqs. I|21|l - (|23|l . without the molec- 
ular field for the experimental^arl^^eters i of a case (ii) 
collapse (acoUapsc = — 12ao), starting with a gaussian ini- 
tial state and no uncondensed atoms due to the absence 
of interactions (ajnit = 0). 

We found that the BEC docs not collapse earlier than 
in the corresponding spherical GPE simulations. The 
number of atoms in excitation modes, shown in Fig. |21 
does not grow fast enough to accelerate the collapse. 
Only just before the collapse time predicted by the GPE 
are large numbers of uncondensed atoms created. We 
compared our generated numbers of uncondensed atoms 
with the predictions of Yurovsky 22] for the central, 
approximately homogeneous region of a large conden- 
sate {Ninit = 50000, Qinit — SOflo). For the case of 
Ocoiiapsc — — 12ao and times much smaller than the non- 
linear interaction time tNL = (^^10(0)1^/^)""'^ ~ 4 ms 
we found that the theories agree. 
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FIG. 2: Number of condensed (a) and uncondensed (b) atoms 
during a case (ii) collapse with acoUapsc ~ — 12ao until t = 
7.3 ms. The initial atom number was No ~ 6000, all in the 
condensate. A three-body loss term with K3 = 1 x 10~^^ 



collapse 



{a/Uo 



-12) is 7.45 ms. Numerical requirements 
prevent the simulations from reaching the exact GPE collapse 
point. Inset of (b): initial evolution of uncondensed-atom 
number. 



Since the main focus of this paper is the investigation 
of ^coiiapso, we did not compute an HFB solution beyond 
the collapse point, where the numerical cost would have 
become immense. The main effort went into the com- 
putation of the correlation functions of the condensate 
fluctuations Gn and Ga from which the uncondensed 
density Gat is derived. Fig. 01 shows the initial evolu- 
tion of the densities of the condensed and uncondensed 
atoms and the corresponding peak densities for the sim- 
ulated period of time. To illustrate the structure of the 
complete correlation functions Gf^/A{'R,~R') some repre- 
sentative samples are displayed in Fig. ^ 

The results presented were computed on a grid with 
Nr = 256 points for R and i?', ranging from -15 ^m 
to -1-15 /im each, and eight Legendre polynomials. We 
checked that the results do not vary for changed spa- 
tial grid size or number of Legendre polynomials. The 
timesteps were 50 ns and 100 ns to check convergence. 
The variation of Nr also verified the cut-off independence 
of our simulations. Our multiprocessor code was run on 
8 or 16 CPUs and employed the RK4IP algorithm 
developed by the BEC theory group of R. Ballagh at the 
University of Otago |^. Parts of the simulations were 
done with the aid of the high level language XMDS [s^l . 



B. Resonance theory 

We also checked that the implementation of the molec- 
ular field as described in section IIIBI (thus an energy- 
dependent model of scattering properties) does not 
change the main result of this paper; that condensate 
depletion is not sufficient to significantly accelerate the 
collapse. If we want to apply the resonance theory for- 
malism to the exact experimental situation, the ini- 
tial state for the molecular field needs to be carefully ad- 
justed. Solving the molecular Eq. H15|) approximately. 
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FIG. 3: a) Time evolution of condensate density n^on^ — \4>\'^ 
and b) density of uncondensed atoms riuucond = Gn for a 
case (ii) collapse with cicollapso ~ —i'2ao- Above we display 
the corresponding peak densities until t — 7.3 ms. 
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FIG. 4: Equal amplitude contours of correlation functions 
Gjv/^(Ro, R) for fixed Ro. We plot correlations for R varying 
in a plane that contains the trap centre and the point Ro. 
The origin of cartesian coordinates (Aa;tang, Axj-adial) in this 
plane is taken to be Ro, a distance |Ro| — Ro from the trap 
centre. The As^adial ^-^is connects Ro with the trap centre, 
which is thus located at the bottom-centre of each plot. The 
contours are equally spaced between the maximum value and 
1/10 (3/10) of it, for Gn (Ga)- Normal density Gn around a) 
Rft — 5.8 fim, b) Ro = 3.2 fim, c) Ro = 0.6 fim. Anomalous 
density Ga around d) Ro — 5.8 fj,m, e) Ro = 3.2 fj,m, f) 
Ro = 0.6 fim. The snapshot is taken during a case (ii) collapse 
with a — — 12ao at revolve ~ 5 nis. 



far away from resonance (a condition fulfilled for the 
collapse-scenarios concerned), yields |17j: 

,^„,(x,t„) = |-,/)2(x,io). (24) 

Simulations done with this particular molecular initial 
state show that the molecular field performs rapid oscil- 
lations around the instantaneous value of g/ (21^)0^ (x, t). 
The effect on the interaction terms in the modified GP 
equation H14|l is just the generation of an effective atom- 
atom interaction of strength Ucg, supporting the view 
presented in that far off resonance the molecular 
field is not important for the condensate evolution. Due 
to the large atom-molecule coupling, g(f>rn is of the or- 
der of J7|0a|^ even for the very small molecular popula- 
tion that we observed. In order to make predictions with 
the same precision as those presented in section Fill Al it 
might therefore be necessary to include quantum fluctu- 
ations of the molecular field as well. Preliminary results 
neglecting these do not indicate production rates of un- 
condensed atoms that differ qualitatively from those in 
Fig. 13 Nonetheless a more careful implementation of the 
molecular field might be a subject for further studies. 



IV. DISCUSSION 

After our investigations an open question remains: 
what is the correct quantitative theory that describes the 
case (ii) collapse experiment? We therefore discuss some 
limitations of our approach, and possible improvements 
and alternatives. 

First it might be possible that a quantum treatment 
with the proper experimental cylindrical symmetry yields 
agreement, but it is not obvious how a changed symmetry 
would cause a qualitative change in the low production 
rates of uncondensed atoms that we found. 

Another factor that might be relevant is any uncon- 
densed population in the initial state, which we assumed 
to be zero. But preliminary simulations that include an 
initial thermal cloud show only a small acceleration of 
the collapse (of the order of 5%) for the experimental 
temperature of 3 nK. Here we estimate the collapse time 
from the onset of rapid decline of the condensate popula- 
tion. For temperatures as high as 7 nK the collapse time 
is about 10% shorter [s^. But even this cannot by it- 
self account for the discrepancy between theoretical and 
experimental data. 
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A drawback of the HFB method is that the interac- 
tion between uncondensed and condensed atoms is not 
properly renormalized |34], since this would re quir e the 
inclusion of higher order correlation functions '3^ l35l| . 
Since the effects of renormalization on the coupling con- 
stant are small in the parameter regime of case (ii), it is 
not a priori clear whether a more complete renormaliza- 
tion would significantly alter the number of uncondensed 
atoms. We applied the gapless HFB method of js^ to a 
case (ii) collapse and found no changes in our results. 

A further step could be to employ a different factorisa- 
tion scheme for the correlation functions, for example the 
method of non-commutative cumulants put forward in 
|3(i| which differs from ours in its truncation method. To 
determine the rate of spontaneous production of uncon- 
densed atoms in that formalism we would need to include 
the higher order correlation functions in our simulations. 

Finally, the growth of the uncondensed fraction by 
elastic collisions as investigated in J^, J^J might also 
play an important role in a case (ii) collapse. Its physics 
seems not to be described by the HFB equations since 
the creation of pairs can be identified as the leading con- 



tribution to the production of uncondensed atoms in the 
limit of weak coupling Uo ■ This does not mean that rele- 
vant physics is missing in the HFB method, as a coherent 
elastic collision process with one atom scattering into a 
previously uncondensed mode might transfer coherence 
into that region of Hilbert space [33. In that case the 
above process would not populate Gn since we define the 
quantum field x to have zero expectation value. 



Acknowledgments 

We would like to thank M. Holland for providing the 
code used in for the purpose of comparison. We also 
thank R. Duine, M. Davis and B. Blakie for discussions, 
as well as S. Metens for his comments. This research 
was supported by the Australian Research Council un- 
der the Centre of Excellence for Quantum- Atom Optics 
and by an award under the Merit Allocation Scheme of 
the National Facility of the Australian Partnership for 
Advanced Computing. 



M. Greiner, O. Mandel, T. Esslinger, T. W. Hansch, and [20 
I. Bloch, Nature 415, 39 (2002). 

E. A. Donley, N. R. Claussen, S. L. Cornish, J. L. 
Roberts, E. A. Cornell, and C. E. Wieman, Nature 412, [21 
295 (2001). 

C. M. Savage, N. P. Robins, and J. J. Hope, Phys. Rev. [22 
A 67, 014304 (2003). [23 
J. L. Roberts, N. R. Claussen, S. L. Cornish, E. A. Don- 
ley, E. A. Cornell, and C. E. Wieman, Phys. Rev. Lett. [24 
86, 4211 (2001). 

J. M. Gerton, D. Strekalov, I. Prodan, and R. G. Hulet, [25 
Nature 408, 692 (2000). 

C. C. Bradley, C. A. Sackett, and R. G. Hulet, Phys. [26 
Rev. Lett. 78, 985 (1997). 
S. K. Adhikari, Physics Letters A 296, 145 (2002). [27 
S. K. Adhikari, J. Phys. B 37, 1185 (2004). 
H. Saito and M. Ueda, Phys. Rev. A 65, 033624 (2002). 
H. Saito and M. Ueda, Phys. Rev. Lett. 86, 1406 (2000). [28 
H. Saito and M. Ueda, Phys. Rev. A 63, 043601 (2001). 
M. Ueda and H. Saito, J. Phys. Soc. Jpn. Suppl. C 72, [29 
127 (2003). 

L. Santos and G. V. Shlyapnikov, Phys. Rev. A 66, [30 
011602(R) (2002). 

W. Bao, D. Jaksch, and P. A. Markowich, J. Phys. B 37, 
329 (2004). [31 
R. A. Duine and H. T. C. Stoof, Phys. Rev. A 68, 013602 
(2003). 

S. Metens, G. Dewel, and P. Borckmans, Phys. Rev. A 
68, 045601 (2003). [32 
E. A. Calzetta and B. L. Hu, Phys. Rev. A 68, 043625 [33 
(2003). [34 
R. A. Duine and H. T. C. Stoof, Phys. Rev. Lett. 86, 
2204 (2001). [35 
J. N. Milstein, C. Menotti, and M. J. Holland, New J. 
Phys. 5, 52 (2003). [36 



N. R. Claussen, S. J. J. M. F. Kokkelmans, S. T. Thomp- 
son, E. A. Donley, E. Hodby, and C. E. Wieman, Phys. 
Rev. A 67, 060701(R) (2003). 

This scenario is described in [T^ but is based on unpub- 
lished data. 

V. A. Yurovsky, Phys. Rev. A 65, 033605 (2002). 

J. K. Chin, J. M. Vogels, and W. Ketterle, Phys. Rev. 

Lett. 90, 160405 (2003). 

J. L. Roberts, N. R. Claussen, S. L. Cornish, and C. E. 
Wieman, Phys. Rev. Lett. 85, 728 (2000). 
J.-P. Blaizot and G. Ripka, Quantum Theory of finite 
systems (MIT Press, 1986). 

C. J. Pethik and H. Smith, Bose-Einstein condensation 
in dilute gases (Cambridge University Press, 2002). 
S. J. J. M. F. Kokkelmans, J. N. Milstein, M. L. Chiofalo, 
R. Walser, and M. J. Holland, Phys. Rev. A 65, 053617 
(2002). 

S. J. J. M. F. Kokkelmans and M. J. Holland, Phys. Rev. 
Lett. 89, 180401 (2002). 

M. Holland, J. Park, and R. Walser, Phys. Rev. Lett. 86, 
1915 (2001). 

Details of the machine are given on the web site of the 
National Faciltity of the Australian Partnership for Ad- 
vanced Computing: http://nf.apac.edu.au/ 
The RK4IP method is described in the Ph.D. the- 
sis of B.M. Caradoc-Davies which is online at: 
http://www.physics.otago.ac.nz/bec2/bmcd/ This is a 
pseudo-spectral method with a Runge-Kutta time step. 
Online at www.xmds.org. 

M. J. Davis and B. Blakie, private communication (2004). 
R. A. Duine and H. T. C. Stoof, J. Opt. B Quant. Semi- 
class. 5, 212 (2003). 

N. P. Proukakis, S. A. Morgan, S. Choi, and K. Burnett, 
Phys. Rev. A 58, 2435 (1998). 

K. Goral, T. Kohler, T. Gasenzer, and K. Burnett, J. 



Mod. Opt. 51, 1731 (2004). 34, 4487 (2001). 

M. J. Davis, R. J. Ballagh, and K. Burnett, J. Phys. B 



